//
//  Move.cpp
//  phi_four
//
//  Created by Snir Gazit on 3/10/12.
//  Copyright (c) 2012 __MyCompanyName__. All rights reserved.
//

#include "Move.h"
// General Move

void Move::init(MTGen* eng, Lattice * lat) {
	eng_ = eng;
	lat_ = lat;
	zeroOne_.init(eng, UniRealParams(0, 1));
}
bool Move::accept(double p) {
	if (p > 1)
		return true;
	double draw = zeroOne_.draw();
	if (p > draw)
		return true;
	else
		return false;
}
bool Move::checkMove() {
	int totalBonds;
	return lat_->testConfiguration(totalBonds);
}
// Shift Move
void ShiftMove::init(MTGen *eng, Lattice *lat, QmcParameters &qmcParameters,
		BoltzWeightType boltzRatio) {
	static_cast<Move *>(this)->init(eng, lat);
	chooseDirection_.init(eng,
			UniIntParams(0, 2 * lat->getCordinationNumber() - 1));
//	createRemove_.init(eng, BinomialVariableParams(1, 0.5));
	//openCloseRatio_ = qmcParameters.openCloseRatio.getVal();
	g_ = 0.5 * qmcParameters.g.getVal();
	auto itR = RemoveRatio_.begin();
	auto itC = CreateRatio_.begin();
	*itR++ = 0.;
	for (auto it = boltzRatio.begin(); it != boltzRatio.end(); it++) {
		*itC++ = ((*it) * g_);
		if ((it) != (boltzRatio.end() - 1))
			*itR++ = (1.0 / ((*it) * g_));
	}
}
void ShiftMove::makeMove() {
	double p;
	// Get Ira
	Lattice::SiteIter ira = lat_->getIra();
	// Choose Direction
	int dir = chooseDirection_.draw();
	// Create or Remove
	//CreateRemove createRemove = (createRemove_.draw() == 0) ? CREATE : REMOVE;
	Site::BondData cBondData = ira->bondArray_[dir];
	if (cBondData.createRemove_ == CREATE) {
		// Get Direction bond Data
		MoveDir moveDir = cBondData.moveDir_;
		// Neighbor Iterator
		int &totalBonds = cBondData.neighbourSite_->totalBonds_;
		int &bondVal = cBondData.bond_->bondVal[moveDir];
		p = CreateRatio_[totalBonds] / double(bondVal + 1);
		if (lat_->isClosed())
			p = p * 2;
		else if (cBondData.neighbourSite_ == lat_->getMasha())
			p = p * 0.5;
		if (this->accept(p)) {
			bondVal++;
			totalBonds++;
			lat_->setIra(cBondData.neighbourSite_);
		}
	} else {
		// Get Direction bond Data
		MoveDir moveDir = cBondData.moveDir_;
		// Neighbor Iterator
		int &bondVal = cBondData.bond_->bondVal[moveDir];
		if (bondVal != 0) {
			int &totalBonds = ira->totalBonds_;
			p = RemoveRatio_[totalBonds] * double(bondVal);
			if (lat_->isClosed())
				p = p * 2;
			else if (cBondData.neighbourSite_ == lat_->getMasha())
				p = p * 0.5;
			if (this->accept(p)) {
				bondVal--;
				totalBonds--;
				lat_->setIra(cBondData.neighbourSite_);
			}
		}
	}
	assert(this->checkMove());
}
// Jump Move

void JumpMove::init(MTGen *eng, Lattice *lat, QmcParameters &qmcParameters_,
		BoltzWeightType boltzRatio) {
	static_cast<Move *>(this)->init(eng, lat);
	//boltzRatio_ = boltzRatio;
	chooseSite_.init(eng, UniIntParams(0, lat->getNumOfSites() - 1));
	boltzRatio_ = boltzRatio;
}

void JumpMove::makeMove() {
	double p;
	Lattice::SiteIter ira = lat_->getIra();
	int cTotalBonds = ira->totalBonds_;
	int jumpSite = chooseSite_.draw();
	Lattice::SiteIter jump = lat_->getSite(jumpSite);
	int jTotalBonds = jump->totalBonds_;
	p = boltzRatio_[jTotalBonds] / boltzRatio_[cTotalBonds - 1];
	if (this->accept(p)) {
		ira->totalBonds_--;
		jump->totalBonds_++;
		lat_->setIra(jump);
		lat_->setMasha(jump);
	}
	assert(this->checkMove());
}

void CreateLoopMove::init(MTGen *eng, Lattice *lat,
		QmcParameters &qmcParameters, BoltzWeightType boltzRatio) {
	static_cast<Move *>(this)->init(eng, lat);
	chooseSite_.init(eng, UniIntParams(0, lat->getNumOfSites() - 1));
	setCoupeling(boltzRatio, qmcParameters.g.getVal());
	asym_ = qmcParameters.asym.getVal();
}
void CreateLoopMove::setCoupeling(BoltzWeightType boltzRatio, double g) {
	g_ = g * 0.5;
	RemoveRatio_[0] = 0.;
	for (int j = 0; j < MAX_BOND; j++) {
		CreateRatio_[j] = boltzRatio[j] * g_;
		if (j != (MAX_BOND - 1))
			RemoveRatio_[j + 1] = 1 / (boltzRatio[j] * g_);
	}
}
int CreateLoopMove::makeMove() {
	int wormLength = 0;
	NCont_.clear();
	currentLoop_.clear();
	//currentLoop_.resize(avWormLength_);
	// Choose First Site
	int jumpSite = chooseSite_.draw();
	Lattice::SiteIter jump = lat_->getSite(jumpSite);
	lat_->setIra(jump);
	lat_->setMasha(jump);
	// Calc Weight of First Site
	getWeights();
	double totalWeight = accumulate(weights_.begin(), weights_.end(), 0.);
	NCont_.push_back(totalWeight);
	// Start Looping
	do {
		wormLength++;
		nextStep();
		getWeights();
	} while (!lat_->isClosed());
	totalWeight = accumulate(weights_.begin(), weights_.end(), 0.);
	NCont_.push_back(totalWeight);
	double pc = NCont_.front() / NCont_.back();
	if (!accept(pc)) {
		eraseLoop();
		eraseWorm_(1);
		//	assert(this->checkMove());
		return 0;
	}
	eraseWorm_(0);
	assert(this->checkMove());
	return wormLength;
}
double calcProb(double lambda, double beta) {
	return pow((lambda / (1. + lambda)), 1. - beta)
			* std::min(1., pow(lambda, beta));
}
void CreateLoopMove::eraseLoop() {
	Lattice::SiteIter cIra = lat_->getIra();
	for (LoopDataIter itl = currentLoop_.begin(); itl != currentLoop_.end();
			itl++) {
		Site::BondData &cBondData = *(*itl);
		Lattice::SiteIter nextIra = cBondData.neighbourSite_;
		// int &totalBonds = nextIra->totalBonds_;
		int &bondVal = cBondData.bond_->bondVal[cBondData.moveDir_];
		if (cBondData.createRemove_ == CREATE) {
			bondVal--;
			nextIra->totalBonds_--;
		} else {
			bondVal++;
			cIra->totalBonds_++;
		}
		cIra = nextIra;
	}
}
void CreateLoopMove::getWeights() {
	auto itW = weights_.begin();
// Choose First Site
	const Site::BondArray &cBonds = lat_->getIra()->bondArray_;
	int cTotalBonds = lat_->getIra()->totalBonds_;
	for (auto it_b = cBonds.begin(); it_b != cBonds.end(); it_b++) {
		auto bond = *it_b;
		int bondVal = bond.bond_->bondVal[bond.moveDir_];
		if (bond.createRemove_ == CREATE) {
			int totalBonds = bond.neighbourSite_->totalBonds_;
			double weight = CreateRatio_[totalBonds] / double(bondVal + 1);
			if (bond.bond_->dir == 2)
				weight *= asym_;
			*itW++ = min(1., weight);
		} else {
			int totalBonds = cTotalBonds;
			double weight = RemoveRatio_[totalBonds] * double(bondVal);
			if (bond.bond_->dir == 2)
				weight /= asym_;
			*itW++ = min(1., weight);
		}
	}
}

int CreateLoopMove::chooseDirection() {
	double totalWeight = accumulate(weights_.begin(), weights_.end(), 0.);
	double draw = totalWeight * zeroOne_.draw(), sum = 0.;
	int chosen = -1;
	for (auto itW = weights_.begin(); sum < draw; itW++) {
		++chosen;
		sum += *itW;
	};
	if (chosen == -1)
		chosen = 0;
	assert(!(chosen < 0));
	assert(!(chosen >= weights_.size()));
	return chosen;
}
void CreateLoopMove::nextStep() {
//	chooseDirection_.setParams(DiscreteVariableParams(weights_));
	int chosenDir = chooseDirection();
	Lattice::SiteIter cIra = lat_->getIra();
//    int spaceDir=chosenDir/(2*DIM);
	Site::BondData &cBondData = cIra->bondArray_[chosenDir];
	Lattice::SiteIter nextIra = cBondData.neighbourSite_;
// int &totalBonds = nextIra->totalBonds_;
	int &bondVal = cBondData.bond_->bondVal[cBondData.moveDir_];
	if (cBondData.createRemove_ == CREATE) {
		bondVal++;
		assert(bondVal < MAX_BOND);
		nextIra->totalBonds_++;
	} else {
		bondVal--;
		assert(bondVal >= 0);
		cIra->totalBonds_--;
	}
	lat_->setIra(nextIra);
	currentLoop_.push_back(&cBondData);
//	assert(this->checkMove());
}
double CreateLoopMove::getEraseProb() {
	return boost::accumulators::mean(eraseWorm_);
}
void CreateLoopMove::setAvWormLength(int avWormLength) {
	avWormLength_ = avWormLength;
}
